%Matlab code to generate figures 4 in the manuscript "Putting the
%significance of spectral peaks on the level: implications for the
%1470-yr peak in Greenland δ18O" by P. Huybers, Feb 2022.  This
%file was run on Matlab release R2020b and requires the signal
%processing toolbox. 

%Examine synthetic scenarios
addpath Toolbox;

clear;
alpha=0.01;
qredo=0;

if qredo,
%%load data
addpath ../Data;
GISP2;
timeo=G(1:end,3)/1000;
Oo=G(:,2);
pl=find(timeo>=20 & timeo<=50);
O=Oo(pl);
time=timeo(pl);
mdt=mean(diff(time));
pl=find(timeo>=0 & timeo<=70);
Oo=Oo(pl);
timeo=timeo(pl);

%interpolate to uniform intervals
ti=time(1):mdt/10:time(end);
Oi=interp1(time,O,ti); 
Oi=smoothPH(Oi,11); 
time=ti(1:10:end); 
O=Oi(1:10:end);
dt=diff(time(1:2));

%%spectral analysis 
opts.dt=dt;
opts.nw=2;
opts.mdlV=1;
dof=2*(2*opts.nw-1);
k=dof/2;
bias=psi(0,k)+log(2)-log(k*2);
conf=log(gaminv(1-alpha,k,1/k));
confb=log(gaminv(1-alpha/(length(time)/dof),k,1/k));
Vexp=psi(1,dof/2);  %expected variance of log(P)-log(Pc)
realizations=10000;
ds=0.1*1/1.470;
so=1/1.470;

%Make realizations
opts.p=1;
opts.q=1;
[dummy s coeffs]=specw(detrend(O),dt,opts.nw,opts);
opts.ar=coeffs.p;
opts.ma=coeffs.q;
mdl = arima('MA',opts.ma,'AR',opts.ar,'Variance',opts.mdlV,'Constant',0);
x=simulate(mdl,length(time),'NumPaths',realizations);

for rep=1:3,
    if rep==1,
        opts.p=1;
        opts.q=0;
        P(rep).opts=opts;
    end;
    if rep==2,
        opts.p=1;
        opts.q=1;
        P(rep).opts=opts;
    end;
    if rep==3,
        opts.p=3;
        opts.q=2;
        P(rep).opts=opts;
    end;

    for iter=1:realizations,
        %null
        xo=detrend(x(:,iter));
        [dummy s]=pmtmPH(xo,opts.dt,opts.nw,0); 
        dummy=log(dummy)-mean(log(dummy))+bias;
        P(rep).init(iter,:)=dummy; 
                
        xpre=prewhitenPH(x(:,iter),opts.p,opts.q);
        xpre=detrend(xpre);
        [dummy s]=pmtmPH(xpre,opts.dt,opts.nw,0); 
        dummy=log(dummy)-mean(log(dummy))+bias;
        P(rep).pre(iter,:)=dummy; 

        %alternate
        do=1.3*sin(2*pi*(time/1.470+rand))';
        xo=xo+do;
        [dummy s]=pmtmPH(xo,opts.dt,opts.nw,0); 
        dummy=log(dummy)-mean(log(dummy))+bias;   
        A(rep).init(iter,:)=dummy; 
        pls=find(s>=so-ds & s<=so+ds);
        fpo(iter,rep,1)=sum(dummy(pls))/sum(dummy);
    
        pl=find(s>0.9/1.470 & s<1.1/1.470);
        dummy=exp(dummy);
        vfrac(rep,iter)=sum(dummy(pl))/sum(dummy);
        
        xpre=prewhitenPH(x(:,iter)+do,opts.p,opts.q);
        xpre=detrend(xpre);
        [dummy s]=pmtmPH(xpre,opts.dt,opts.nw,0); 
        dummy=log(dummy)-mean(log(dummy))+bias;
        A(rep).pre(iter,:)=dummy; 

        x1=P(rep).pre(:);
        nx1=mean(x1>conf);
        x2=max(P(rep).pre,[],2);
        nx2=mean(x2>confb);

        x1=A(rep).pre(:);
        ax1=mean(x1>conf);
        x2=max(A(rep).pre,[],2);
        ax2=mean(x2>confb);

        [rep iter nx1 nx2 ax1 ax2 vfrac(rep,iter)],    
    end;
end;
save GISPsynth.mat;
else,
    load GISPsynth.mat x P A s k time opts Vexp;
end;

%%Statistics used for plotting
alpha=0.01;
conf=log(gaminv(1-alpha,k,1/k));
confb=log(gaminv(1-alpha*2*k/length(time),k,1/k));

so=1/1.470;
[dummy pls]=min(abs(s-so)); pls=pls+2;  %Find 1/1470 year peak
for rep=1:3,
    x=P(rep).pre;
    v=var(x'); 
    pl=find(v>=inf);
    v(pl)=NaN;
    x(pl,:)=NaN;
    nv1(rep)=nanmean(v);
    nF1(rep)=nanmean(v/Vexp);
    nF1_001(rep)=prctile(v/Vexp,1);
    nF1_99(rep)=prctile(v/Vexp,99);

    v=var(exp(x'));
    nF1b(rep)=nanmean(v);

    nx1(rep)=nanmean(x(:,pls)>conf); %H0 coverarge
    x2=max(x,[],2);
    nx2(rep)=nanmean(x2>confb); %H0max coverarge

    x=A(rep).pre;
    v=var(x'); 
    av1(rep)=nanmean(v);
    aF1(rep)=nanmean(v./Vexp);

    v=var(exp(x'));
    aF1b(rep)=nanmean(v);

    ax1(rep)=nanmean(x(:)>conf); %H1 coverarge
    x2=max(x,[],2);
    ax2(rep)=nanmean(x2>confb); %H1max coverarge
end;

%%Figure including null and alternate distributions
figure(4); clf; 
[dummy pls]=min(abs(s-so)); pls=pls+2;
load EX_GISP.mat mag1470 vP;
xi=-1.5:0.01:3;
conf=log(gaminv((1-alpha),k,1/k));
confb=log(gaminv(1-alpha*2*k/length(time),k,1/k));
for rep=1:3,
    if rep==1, 
        txt={'(a)'; ' '; 'ARMA(1,1)'; 'ARMA(1,0)'};
        txt2='(b)';
        Pobs=mag1470(1);
        col=[1 0.1 0.1];
    end;
    if rep==2, 
        txt={'(c)'; ' '; 'ARMA(1,1)'; 'ARMA(1,1)'};
        txt2='(d)';
        Pobs=mag1470(2);
        col=[0.1 0.1 1];
    end;
    if rep==3, 
        txt={'(e)'; ' '; 'ARMA(1,1)'; 'ARMA(3,2)'};
        txt2='(f)';
        Pobs=mag1470(3);
        col=[0.1 0.7 0.5];
    end;
    subplot(3,4,[1 2 3]+(rep-1)*4); cla; hold on; 
    yv=var(P(rep).pre,[],2);
    plv=find(yv<=inf); 
    vpercent(rep)=mean(yv<=Vexp);
    
    y=P(rep).pre(:,pls);
    xc(1,rep)=xcPH(y,yv,1);
    H0=ksdensity(y(plv),xi); 
    Palpha=prctile(y,(1-alpha)*100);
    
    y=max(P(rep).pre,[],2);
    H0b=ksdensity(y(plv),xi);
    xc(2,rep)=xcPH(y(plv),yv(plv),1);
    Palphab=prctile(y,(1-alpha)*100);

    y=max(A(rep).pre,[],2);
    H1=ksdensity(y(plv),xi); 
    xc(3,rep)=xcPH(y(plv),yv(plv),1);

    plot(xi,H0,'k--','linewidth',1.5); 
    plot(xi,H0b,'k','linewidth',1.5); 
    plot(xi,H1,'color',[0.4922    0.6211    0.6484],'linewidth',1.75); 
    axis tight; h4=axis;
    plot(confb*[1 1],[0 h4(4)],'color',[0.8672    0.4531    0.0859],'linewidth',1.5);
    plot(Palphab*[1 1],[0 h4(4)],'k','linewidth',1.5);
    plot(Pobs*[1 1],[0 h4(4)],'color',col,'linewidth',1.5);
    font(gca,13);
    h=xlabel('log spectral density','fontsize',13);
    h=ylabel('normalized histogram','fontsize',13);
    h=text(-1.4,h4(4)*0.7,txt,'fontsize',12);

    subplot(3,4,rep*4); cla; hold on; 
    vi=0:0.01:3;
    ks=ksdensity(yv/Vexp,vi); 
    plot(vi,ks,'k','linewidth',1.5);
    axis tight; h4=axis;
    plot([1 1],h4(3:4),'color',[0.8672    0.4531    0.0859],'linewidth',1.5); 
    plot(vP(rep)/Vexp*[1 1],h4(3:4),'color',col,'linewidth',1.5); 
    h=xlabel('variance ratio, F','fontsize',13);
    h=ylabel('normalized histogram','fontsize',13);
    h=text(0.1,h4(4)*0.875,txt2,'fontsize',12);

end;

return;
print -depsc fig4_revised.eps;
!pstopdf fig4_revised.eps
!open fig4_revised.pdf


